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COLLECTIVE MOTIONS OF A ONE -DIMENSIONAL 
SELF-GRAVITATING SYSTEM 

By Frank Hohl and Janet W. Campbell 
Langley Research Center 

SUMMARY 

A one -dimensional computer model for a collisionless, self- gravitating system is 
used to investigate the mixing phase during which a stellar system approaches a quasi- 
equilibrium state. The results of the calculations are compared with the theory of the 
mixing phase proposed by Lynden-Bell. Previous investigations have shown that for 
simple (uniform phase density) initial conditions the quasi-equilibrium state, which is 
reached after several crossing times, is close to the distribution predicted by Lynden- 
Bell. In the present paper more complicated initial conditions (two-phase density sys- 
tems) are used to test the Lynden-Bell theory. For most systems investigated, the final 
state is close to the Lynden-Bell distribution. One exception is noted. 

INTRODUCTION 

In considering the evolution of stellar systems, the two fundamental time scales of 
interest are the crossing time t c and the thermalization time T r . The crossing time 
is associated with the collective motions of the stellar system and the thermalization 
time is associated with the binary interactions (collisions) of the individual stars. The 
crossing period is 

T - 

C \/47tG p 

where G is the gravitational constant and p is the mass density. The mixing phase 
during which a collisionless system tends to approach a quasi-equilibrium state is of the 
order of r c . For the systems considered in the present report, the thermalization time 
is much larger than the crossing time, r r » r c . Thus, the mixing phase is well sepa- 
rated from the thermalization phenomena. 

Recently, Lynden-Bell (ref. 1) proposed a theory of the mixing phase during which a 
collisionless stellar system approaches a quasi-equilibrium state. By using a new type 
of statistics related to Fermi-Dirac statistics, he derived a final or quasi -equilibrium 
distribution. In the present paper, computer experiments with a one-dimensional sheet 



model (refs. 2 to 4) are performed to determine how closely a collisionless self- 
gravitating system will approach the Lynden-Bell distribution. The case for the most 
simple initial conditions (uniform phase density) has previously been investigated (ref. 5). 
These results showed that from 80 percent to 96 percent of the stars in the system 
approached the Lynden-Bell distribution. Similar results for the uniform initial phase 
density case were obtained by Cohen and Lecar (ref. 6). Henon (ref. 7) used a model of 
concentric spherical mass shells to test the Lynden-Bell theory for the simple uniform- 
phase -density case and also obtained essentially the same results. It was found previ- 
ously (ref. 8) that for the case of a plasma, an initially uniform phase -density system 
will approach the Lynden-Bell distribution. 

For the simple case of an initially uniform phase -density distribution, the final 
state of the system is close to the distribution predicted by Lynden-Bell. It is of interest 
to determine whether more complicated initial distributions also approach the Lynden- 
Bell distribution. 

The present investigation includes initial conditions having two different phase 
densities. Also, effects of the presence of different mass groups are investigated. 

SYMBOLS 


E gravitational field 

f(x,v) distribution function, phase space density 

F(e) energy distribution function 

G gravitational constant 

m mass or mass per unit area 

N total number of stars in system 

P potential energy 

T kinetic energy 

t time 

v velocity 


2 



X 


position coordinate 


/3 1//3 corresponds to a "temperature" 

e total energy per unit mass 

r) initial phase density 

p. "Fermi energy" 

t c crossing time 

t c q crossing time at t = 0 

T r thermalization time 

cp gravitational potential 

< > averaged quantity 

Subscripts: 

i,j summation indices 

1,2 describes two different phase -density regions 

STATISTICS 


The equations describing a one -dimensional collisionless system are the collision- 
less Boltzmann equation 


9f 9f 

— + v — 
9t dx 


-^ 9 A = 0 

dx 9V 


( 1 ) 


and the Poisson equation 


v2 £ p = 47rG ^ f dv (2) 

The most probable final state for a system described by equations (1) and (2) was recently 
discussed by Lynden-Bell (ref. 1). He assumes the phase space to be divided into a large 
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number of elements of phase fp He then argues that because of the large number of 
equivalent elements and the violent changes in the mean field, any distinct distribution 
of elements of phase is equally likely to occur subject to the conditions that the total 
energy and the number of elements of phase with a given initial density are conserved. 

The elements of phase fj[ are distinguishable. According to equation (1), the elements 
of phase cannot overlap and, therefore, follow an exclusion principle. Such a phase fluid 
should obey a fourth type of statistics which, in general, is different from that of Maxwell- 
Boltzmann, Bose -Einstein, and Fermi-Dirae. By evaluating the most probable distribu- 
tion of elements of phase, Lynden-Bell obtains the coarse-grained (macroscopic) distri- 
bution 


j 1 + 


exp Oi( 6 - pjj] 


I exp[-^(e - pj) 


(3) 


} 


where the summation occurs because initially there are elements of phase with different 
densities 77 j . 

For the special case where the initial f has a constant value rj over certain 
regions of phase space and is zero outside these regions (water-bag distribution), the 
distribution reduces to 


<f> 


V 

1 + exp[/3(e - p)] 


(4) 


where' e is the total energy of a star and the two constants /3 and p are determined 
by the conservation of energy and area in phase space. The distribution given by equa- 
tion (4) is formally identical with the Fermi-Dirac distribution. The parameter p 
corresponds to the Fermi energy, whereas the parameter 1//3 corresponds to a 
"temperature" and is a measure of the excitation above the minimum energy state corre- 
sponding to 1//3 = 0. This minimum-energy state is discussed by Hohl and Feix (refs. 4 
and 9) for the case of a one-dimensional system and in this paper is referred to as the 
stationary state. 

Equation (4) can be written in the form 


/3(e - p) = log 


(aii\ 

U-<f> 


(5) 


so that a plot of log e |( f )j [r] - ( f )) as a function of energy should give a straight line. 
For a highly nondegenerate system such as that investigated by Buneman (ref. 2) ( f ) « 77 , 
and the distribution given by equation (4) will approach a Maxwellian. On the other hand, 
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a distribution near the stationary state of a single contour water -bag distribution (refs. 4 
to 6) corresponds to a low-temperature (1//3 ~ 0) Fermi distribution. For the case of a 
system with two different initial-phase densities tjj_ and 7? 2 , the following relationships 
hold (appendix I of ref. 1) 


< f l > = 


exp[-/3i(e - Ml) 


1 + exp 

'-Jl\ 

<e - Ml)J + exp[-/3 2 | 

[ e - M2)_ 


( 6 ) 


and 


< f 2 > = V2 


exp[-/3 2 (e - M 2 )] 


1 + 


exp[]/3 1 (e - mi)] + exp[-/3 2 (e - M 2 ). 


(7) 


or 


and 


h{ e 


Ml) = log e 



(fl) 
<£i > 

Vi 



h( e 


m 2 ) = log e 



<f2) 

(fl) ( f 2 >\ 
VI V 2 j 


( 8 ) 


(9) 


THE MODEL 

A one -dimensional sheet model (refs. 2 to 4) was used for the computer simulation 
of the mixing phase of collisionless systems. Consider a system of N mass sheets with 
ordered sheet positions xi,x 2 ,X 3 ,. . .,x n such that xj=Xj + i. The gravitational field Ej 
acting on a sheet at Xj is obtained from the net mass to the left of Xj; that is, 

N 

Ej = 2 ttG I mi sgn(xj - x^ (10) 

i=l 

where is the mass per unit area of the sheet at x^ and 

sng(x) = -1 
sgn(x) = 0 
sgn(x) = 1 
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(x > 0) 

(x = 0) > (11) 

(x < 0) 

J 



The velocity and position of the sheets are advanced stepwise in time by using the time- 
centered finite -difference equations 


( 12 ) 

(13) 

The process of ordering the sheets, calculating the field, and advancing velocities and 
position by a small time step 6t is repeated until the desired evolution of the system is 
achieved. 

RESULTS AND DISCUSSION 



The Lynden-Bell theory of violent relaxation (ref. 1) does not account for any 
effects due to the presence of different mass groups. It is therefore important to ensure 
that no mass segregation occurs when different mass groups are present. The effect of 
the presence of two different mass groups was determined by repeating some of the calcu- 
lations presented in reference 5, but by representing the upper stream (positive velocity 
stars) by 2000 mass sheets each of mass 1 (per unit area) and the lower stream by 
1000 mass sheets each of mass 2. The resulting evolution of the system in phase space 
is shown in figure 1. The time is shown in crossing periods t c q calculated for the 
initial density. The initial distribution is equivalent to that of the system shown in fig- 
ure 1 of reference 5. The evolution of the system shown in figure 5 of reference 5 was 
also repeated by replacing the upper stream by 2000 stars of mass 1/2 and the lower 
stream by 1000 stars of mass 1. Figure 2 shows the resulting evolution in phase space. 
The evolution displayed in figures 1 and 2 is practically identical to that shown in fig- 
ures 1 and 5 of reference 5. Also, when log e [(f) /(rj - (f)jj is plotted separately for 
the light and heavy stars, a practically identical distribution is obtained for both cases. 
Each one of the plots of loggj^f)/^ -(f)) gives the same variations as the corre- 
sponding plot presented in reference 5 for the same initial phase-density system. Thus, 
the final distribution when two different mass groups are present is the same as that for 
the case when all sheets are of equal mass. 

In figure 3 the percent kinetic-energy difference given by 


i ^ m jV^ (light stars) - irqvj^(heavy stars) 


~ ^ mp/^ (light stars) + mjVj2(heavy stars) 
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is plotted as a function of time for the two systems shown in figures 1 and 2. It can be 
seen that in spite of the large fluctuations, there is no apparent systematic increase in 
the kinetic energy of the light stars; that is, no equipartition of energy occurs on a time 
scale r c . These results are to be expected since Hohl and Broaddus (ref. 10) have found 
that the thermalization or relaxation effects for the one -dimensional model occur on a 
time scale 

r r ~ N 2 t c (15) 

where N is the number of mass sheets in the system. Thus, during the process of 
violent relaxation as simulated by the one -dimensional model, purely collective inter- 
actions of the stars cause the system to approach its quasi -equilibrium state. 

The evolution in phase space of the first two-phase-density systems investigated is 
shown in figure 4. Initially, the central region of the system bounded by x = ±187.5 and 
v = ±375 contains 1008 stars of mass 1. The two outer regions bounded by x = ±187.5, 
x = ±375, and v = ±375 each contain 1008 stars of mass 1. Therefore, the initial den- 
sity in the central region defined by is one-half that of the two outer regions defined 

by f 2 - 

To compare the final state of the system with the distribution given by equations (6) 
and (7), the phase space is divided into a number of cells. Each mass sheet represents a 
point in a cell in phase space. The average number of mass sheets per cell and the 
corresponding energy e = ^-mv2 + mcp are determined by averaging over several dis- 
tributions near the final time t = 47.7 t c q as shown in figure 4. These calculations are 
done separately for stars corresponding to fi and to f 2 - The cell size chosen was 
such that initially, the phase densities (mass per cell) corresponding to fj and f 2 
were pj = 20.16 and P 2 = 40.32. 

Figure 5 shows the variation of the distributions (fi), (f 2 ), and (f)=(fj_ + f 2 ) 
as a function of energy. Also shown as inserts in the same figure are Fj_, F 2 , and F 
which are the distributions of the mass per energy interval corresponding to (fj), (f2)> 
and (f), respectively. The normalizations 

J f dx dv = 1 
^ fl dx dv = 1 
y f 2 dx dv = 1 
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and 


J F de = 3 X 10 4 
j* Fi de = 3 X 10 4 
y F 2 de = 3 X 10 4 


are used throughout the present paper. The dots shown in figure 5 are the numerical 
results obtained from the computer calculations. The solid curves represent the Lynden- 
Bell distribution which gives the best fit to the numerical data; that is, the quantity 


fi (theory) - f^ (numerical) 


is minimized. 


For high energies the value of (f) 


is 


i 

too large and in obtaining the fit the high-energy bump in the distribution was neglected. 
This high-energy bump in the distribution was also obtained by Henon in his spherical 
shell model (ref. 7). Lynden-Bell (ref. 11) suggests that this deviation occurs because 
the high-energy stars have a large period and remain outside the main system for a 
longer time than the low-energy stars. When the high-energy stars return to the center 
of the system, the system has nearly reached a steady state and the mechanism for phase 
mixing would no longer be effective. The values of the constants j3j, fy, Ml, and ^2 
corresponding to the "best" fit shown in figure 5 are also given in figure 5. One of the 
difficulties of the Lynden-Bell distribution is that the parameters of the final state cannot 
be calculated in advance. Nevertheless, the results presented in figure 5 show that the 
final state can be closely approximated by the Lynden-Bell distribution. The dip in (f2> 
obtained from the numerical results can be explained in terms of the initial distribution. 
Initially, the central part has a low phase density and as the system evolves, the outer 
higher phase density tries to displace the lower central density. Since the two-phase 
densities cannot pass through each other (exclusion principle), some portions of the 
lighter region are trapped near the center of the system, and thereby stars corresponding 
to the heavier phase density are prevented from occupying this region. The result is a 
lower phase density in the central region as shown in figure 5. In the previous investiga- 
tion (ref. 5) of single-phase density systems, equation (5) was used for comparing the 
computer results with the Lynden-Bell distribution. An attempt was made to use equa- 
tions (8) and (9) for the present two-phase density cases. However, the numerically 
obtained deviations from the Lynden-Bell distribution for small energies (as shown in 
fig. 5) caused large variations in the logarithmic terms and the results were misleading. 
For this reason equations (6) and (7) are used. 
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The evolution of the kinetic energy for the system shown in figure 4 is given in 
figure 6. Initially, the total potential energy P is nearly 10 times the total kinetic 
energy T. However, after a short time, the virial theorem is satisfied (if the average 
is taken over the oscillation in the energies) and 2T = P. 

The evolution in phase space for the next case is shown in figure 7. Initially, the 
region bounded by x = 0, x = -250, and v = ±312.5 contains 1000 stars of mass 1.6, and 
the region bounded by x = 0, x = 250, and v = ±312.5 contains 1000 stars of mass 0.4. 
Thus, the left-hand region corresponding to f2 has an initial phase density four times 
as large as the right-hand region which corresponds to fj. Also, the heavy stars are 
represented by small rectangles and the light stars by small circles. It can be seen that 
most of the high-energy stars are light stars. The corresponding distribution functions 
obtained near t = 47.7t c q are shown in figure 8. For this case the "best" fit was 
obtained by including all the high-energy bumps. The total distribution (f) and also 
(f 2 ) are fairly well represented by the Lynden-Bell distribution. However, the two 
peaks in (fl) are probably caused again by a combination of trapping between regions 
of (fl) and by nonmixing of high-energy stars. The evolution of the kinetic energy for 
the same system is shown in figure 9. Initially, the potential energy is four times the 
kinetic energy but the energies quickly reach the point where they oscillate around 
P = 2T. 

The next case investigated is one where a region of heavy phase density completely 
encloses a region of low phase density. The evolution in phase space for this system is 
shown in figure 10. The central rectangle defined by x = ±177 and v = ±177 contains 
1000 stars each of mass 0.5. The outer region exterior to the central rectangle and 
enclosed by x = ±250 and v = ±250 contains 1000 stars each of mass 1.5. The two 
regions have equal areas so that the outer region corresponding to fg has a phase den- 
sity which is three times as large as the phase density in the central region which corre- 
sponds to fi. Figure 10 shows that nearly all the high-energy stars are heavy stars 
represented by small rectangles (corresponding to f2). This result was to be expected 
since stars from the inner regions (corresponding to fij cannot penetrate the outer 
region f2- Thus, again the light region is trapped near the center of the system. 

Another point to note is illustrated by the stage of the evolution at t = 3.2 t c q in fig- 
ure 10. It can be seen that parts of the outer region have condensed into two clusters 
which displace some of the lighter regions away from the central part of the system. 
However, these displaced regions are still trapped inside regions of f 2 - Figure 11 
shows a comparison of the distribution functions of this system with the Lynden-Bell 
distribution after the system has reached the quasi -equilibrium state. The results are 
essentially the same as those for the previous two cases; that is, most of the stars follow 
the Lynden-Bell distribution. The variation of the kinetic energy for the system is shown 
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in figure 12. Initially, the potential energy of the system is almost four times as large 
as the total kinetic energy. 

The final system investigated is shown in figure 13. The top and the bottom rec- 
tangle are defined by x = ±83.3, v = ±83.3, and v = ±250 and each of the two rectangles 
contains 506 stars of mass 0.4 represented by a circle. The right and left rectangles 
are defined by x = ±83.3, x = ±250, and v = ±83.3 and each of the two rectangles con- 
tains 506 stars of mass 1.6 represented by a small rectangle. The region of phase space 
covered by the right and left rectangle corresponds to f 2 and has a phase density four 
times as large as the top and bottom rectangles which correspond to f j. Compared with 
the previous three cases, the system has a rather large initial ratio of potential to kinetic 
energy; the initial potential energy is 23 times as large as the total initial kinetic energy. 
Figure 14 shows the final distribution functions for the system obtained near t = 47.7t c q. 
The variation of the kinetic energy for this system is shown in figure 15. For this sys- 
tem it is not possible to obtain any reasonable fit of the numerical results with the Lynden- 
Bell distribution. The relaxation is certainly very violent as required by the Lynden-Bell 
theory. 

However, from the dynamics of the system as displayed in figure 13 it can be seen 
that the final state of the system is rather purely mixed. Initially the two heavy-phase 
density rectangles tend to fall toward the center of the system and most of the stars 
corresponding to the low-phase-density regions are displayed toward higher energies. 

In fact, up to a time t = 4.0t c q, the system rotates (in phase space) like a binary system 
with the heavy-phase -density regions continually approaching each other. At the same 
time, the low -phase -density regions are being pushed to higher energies. At t = 8.0t c q, 
the two heavy regions have merged and occupy the central region of the system. Only 
very few of the stars corresponding to the initial-low-phase-density regions are near the 
center of the system. This state of the system corresponds, of course, to the distribu- 
tions shown in figure 14. These results indicate that if the final state is to be described 
by the Lynden-Bell distribution, the evolution of the system must proceed with sufficient 
interpenetration of regions of different phase density. 

For the single-phase-density case, the stationary state is always stable and is 
described by a zero-temperature (1//3 ~ 0) Fermi-Dirac distribution. In studying the 
evolution of such single-phase-density systems, it was found (refs. 4 and 9) that within 
the limitations of energy conservation, the system approaches the stationary state. Since 
the Lynden-Bell distribution reduces to the Fermi-Dirac distribution for this simple case, 
one would expect the final state of the system to be approximated by the Lynden-Bell dis- 
tribution. Consider now the case of a multiphase density system (refs. 9 and 12). A 
stable stationary two-phase-density system is described by 
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fl(e) = Ml 


(0 e g jx) 


fl(e) = 0 


(e > n) 


and 


* 2 (e) =V2 (Ml < e ^ M2) 

f 2 (e) =0 (e < Mi; e > M 2 ) 


where Ml > M2- & 1 terms of Fermi-Dirac distribution, and f2 can be written as 


f l( e> . m exg [:^ - ^)L. 
1 1 + exp[-/3 1 (e - MlJ 


(16) 


and 


f 2 (e) = ^2 


f exp [- 

■£ 2 ( e “ m 2 )] ex p[- 

-£l(‘ 

; - 

|l + exp 

)[-/3 2 (e - M 2 J ' 1 + exp 

>[-/3; 

i( 6 


(17) 


where l//3i ~ 0 and l/j^ ~ 6- Equations (16) and (17) give a better approximation of 
the final state of the system shown in figures 13 and 14 than the Lynden-Bell distribution. 


CONCLUDING REMARKS 


Since the motion of "stars" in the two-dimensional phase space is always bounded, 
the collisionless, one-dimensional model is well suited for a check of the Lynden-Bell 
theory. It was found that no mass segregation occurs when two different mass groups 
are present. 

Most of the calculations with two different initial phase densities give results in 
good agreement with the Lynden-Bell distribution. The small deviations in these cases 
can be explained by considering the dynamics of the stars. For one case investigated, 
however, the final state of the two phase densities cannot even be approximated by the 
Lynden-Bell distribution. Since the Lynden-Bell distribution is obtained without actually 
considering the dynamics of the stars, it is surprising that for most of the systems 
investigated, the final state does correspond to a Lynden-Bell distribution. One of the 
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disadvantages of the Lynden-Bell distribution is that the parameters of the final state 
cannot be calculated in advance. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., September 19, 1969. 
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